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We give a fairly comprehensive review of wavelets and of their application to 
density-functional theory (DFT) and to our recent application of a wavelet-based ver- 
sion of linear-response time-dependent DFT (LR-TD-DFT). Our intended audience is 
quantum chemists and theoretical solid-state and chemical physicists. Wavelets are a 
Fourier-transform-like approach which developed primarily in the latter half of the last 
century and which was rapidly adapted by engineers in the 1990s because of its advan- 
tages compared to standard Fourier transform techniques for multiresolution problems 
with complicated boundary conditions. High performance computing wavelet codes 
now also exist for DFT applications in quantum chemistry and solid-state physics, no- 
tably the BigDFT code described in this chapter. After briefly describing the basic 
equations of DFT and LR-TD-DFT, we discuss how they are solved in BigDFT and 
present new results on the small test molecule carbon monoxide to show how BigDFT 
results compare against those obtained with the quantum chemistry gaussian-type or- 
bital (GTO) based code DEMon2k. In general, the two programs give essentially the 
same orbital energies, but the wavelet basis of BigDFT converges to the basis set limit 
much more rapidly than does the GTO basis set of DEMon2k. Wavelet-based LR-TD- 
DFT is still in its infancy, but our calculations confirm the feasibility of implementing 
LR-TD-DFT in a wavelet-based code. 
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1. Introduction 

The broad meaning of "adaptivity" is the capacity to make something work better by al- 
ternation, modification, or remodeling. Concepts of adaptivity have found widespread use 
in quantum chemistry, ranging from the construction of Gaussian-type orbital (GTO) ba- 
sis sets, see e.g., the development of correlation consistent bases H] 121 [3, to linear scaling 
methods in density functional theory (DFT) |@J|5]|6l|7]|8]|9l, selective configuration interac- 
tion (CI) methods iflOl [TTJ and local correlation methods based on many-body perturbation 
theory or coupled cluster (CC) theory lfT2l[L3l . This chapter is about a specific adaptive tool, 
namely wavelets as an adaptive basis set for DFT calculations which can be automatically 
placed when and where needed to handle multiresolution problems with difficult boundary 
conditions. 

Let us take a moment to contrast the wavelet concept of adaptivity with other types of 
adaptivity. In other contexts, the adaptive procedure is typically based on a combination of 
physical insights together with empirical evidence from numerical simulations. A rigorous 
mathematical justification is usually missing. This may not be surprising: Familiar con- 
cepts lose a lot of their original power if one tries to put them in a rigorous mathematical 
framework. Therefore, we will not shoulder the monumental and perhaps questionable task 
of providing a rigorous mathematical analysis of all the adaptive approaches used nowadays 
throughout quantum chemistry. Instead we will concentrate on the mathematical analysis 
of a particular electronic structure method which lends itself to a rigorous mathematical 
analysis and application of adaptivity. In contrast with other adaptive methods, multires- 
olution analysis (MRA) with wavelets can be regarded as an additive subspace correction 
and their wavelet representations have a naturally built-in adaptivity which comes through 
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their ability to express directly and separate components of the desirable functions living 
on different scales. 

This combined with the fact that many operators and their inverses have nearly sparse 
representations in wavelet coordinates may eventually lead to very efficient schemes that 
rely on the following principle: Keep the computational work proportional to the number 
of significant coefficients in the wavelet expansions of the searched solution. As there are 
a lot of different wavelet bases with different properties (length of support, number of van- 
ishing moments, symmetry, etc.) in each concrete case we can choose the basis that is most 
appropriate for the intrinsic complexity of the sought-after solution. This fact makes the 
wavelet-based schemes a very sophisticated and powerful tool for compact representations 
of rather complicated functions. The expected success of wavelet transforms for solving 
electronic structure problems in quantum mechanics are due to three important properties: 

(a) the ability to choose a basis set providing good resolution where it is needed, in those 
cases where the potential energy varies rapidly in some regions of space, and less in others; 

(b) economical matrix calculations due to their sparse and banded nature; and (c) the abil- 
ity to use orthonormal wavelets, thus simplifying the eigenvalue problem. Of course, this 
might lead to adaptive methods which are fully competitive from a practical point of view, 
for example, working with a systematic basis instead of GTO bases requires from the on- 
set larger basis sets and the benefit of systematic improvement might be a distant prospect. 
However, we have the more realistic prospect that our rigorous analysis provides new and 
hopefully enlightening perspectives on standard adaptive methods, which we reckon cannot 
be obtained in another way. 

On the otherhand advances in computational technology opened up new opportunities 
in quantum mechanical calculation of various electronic structures, like molecules, crystals, 
surfaces, mesoscopic systems, etc. The calculations can only be carried out either for very 
limited systems or with restricted models, because of their great demand of computational 
and data storage resources. Independent particle approximations, like the Hartree-Fock 
based |[T4l [151 [T6l [171 algorithms with single determinant wave functions, leave out the 
electron correlation and need operation and storage capacity of order iV 4 , if N is the total 
number of electrons in the system. If inclusion of the electron correlation is necessary, CI or 
CC methods can be applied, with very high demand of computational resources (0(iV 6 ) to 
O(Nl)). An alternative way is to use MBPT. The second order perturbation calculations can 
be carried out within quite reasonable time and resource limits, but the results are usually 
unsatisfactory, they just show the tendencies, while the 4th order MBPT needs 0(N 7 ) 
to 0(N 8 ) operations. All these algorithms use the iV-electron wave function as a basic 
quantity. 

Another branch of methods use electron density as the primary entity. Pioneers of this 
trend, like Thomas lfl8l . Fermi lfT9ll20ll . Frenkel ETI and Sommerfeld [22] developed the 
statistical theory of atoms and the local density approximation (LDA). The space around 
the nuclei is separated into small regions, where the atomic potential is approximated as 
a constant, and the electrons are modeled as a free electron gas of Fermi-Dirac statistics 
ll23l l24l |2~T1. Dirac included electron correlation 11251 . which improved the results. Af- 
ter the Hohenberg-Kohn theorems had appeared |[26l . and Kohn and Sham had offered a 
practically applicable method ll27l based on their work, many scientists were motivated to 
work on the theory, and DFT developed into one of the most powerful electronic structure 
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methods. 

Despite the success of density functional theory, it has some drawbacks. The exact 
formula of the exchange-correlation potential is not known, thus chemical intuition and 
measured data are necessary in order to approximate it, and the kinetic energy functional is 
hard to calculate. Powerful approximating formulas are available (see, e.g. Il28l0 . like the 
Thomas-Fermi functional based gradient and generalized gradient expansions, where the 
energy functionals are expressed as a power series of the gradient of the density (the first 
such suggestion was ll29l .) 

Considering the historical development of sophisticated iV-electron methods, a typical 
trend can be observed. Starting with a very simple model, new details are introduced in 
order to improve the results. This scheme is followed in the linear muffin tin orbital method 
(LMTO) QUI where the interatomic regions is replaced by the spherical orbital of an atomic 
potential around the nuclei. Similarly, the linearized augmented plane wave method (APW) 
lOTI and the plane wave pseudopotential approach P2l describe the details of the crystal 
potential differently in different spatial regions. Although they are rather successful, for 
applying any of these models, chemical intuition is needed, free parameters, like the radius 
of the bordering sphere between the two types of potentials, and the boundary conditions 
have to be set. A systematic method, which can handle the different behaviors of the elec- 
tron structures at different spatial domains, or either at different length scale t33l . is the 
longterm requirement of any physical chemists. 

Multiresolution or wavelet analysis, this rapidly developing branch of the applied math- 
ematics, is exactly the tool for statisfy all the need of any chemical physicists/physcial 
chemists. From mathematical point of view wavelet analysis is a theory of a special kind 
of Hilbert space basis sets. Basis sets are commonly used in all electron structure calcu- 
lations, as the wave function is usually expanded as linear combination of some kind of 
basis functions. Thus the operator eigenvalue problem is reduced to an algebraic matrix 
eigenvector problem. The resulting algebraic equations are easier to solve, well known 
algorithms and subroutine libraries are available, however, the difficulty of choosing the 
proper basis set arises. If linear combination of atomic orbitals (LCAO) is used, the atomic 
basis functions are Slater or Gaussian-type of functions P4l l35l . the selection of atomic 
orbitals needs chemical intuition, which is a result of long time's experience, and can not be 
algorithmized. Both basis sets are non-orthogonal, and lack the explicit convergence prop- 
erties ll36l . Moreover, calculation of operator matrix elements with Slater-type orbitals is 
complicated, their integrals have to be treated numerically. Although integrals of Gaussian 
functions are analytically known, the Gaussisn-type basis does not reflect the nuclear cusp 
condition of Kato 071 . which reflects on singularities of the iV-electron wave function in 
the presence of Coulomb-like potentials. Since then it turned out that for high precision nu- 
merical calculations it is essential to satisfy these requirements. However, while the nuclear 
cusp condition is relatively easy to fulfill by Slater-type orbitals (STO), the electron-electron 
cusp is extremely hard to represent. In general, GTO-based/STO-based DFT codes gives 
reliable results with a relatively small number of basis functions, making them optimal for 
large scale computations where high accuracy is less crucial. On the other hand there is no 
consistent way to extend these basis sets and thereby converge the results with respect to 
the size of the basis. 

The second type of basis set covers the system-independent functions such as plane 
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waves ll32l or wavelets ll38l . The main advantable of these basis sets, is that their size can 
be systematically increased until the result of the calculation has converged, and are gen- 
erally considered to be more accurate than the former type. The number of basis function 
required to obtain convergence is normally so large that direct solution of the matrix eigen- 
value problem within the entire basis space is not possible. Instead one has to use iterative 
methods to determine the lowest (occpied) part of the spectrum ll32l . In solid state physics, 
where more or less periodic systems are studied, choosing plane wave basis sets is rather 
usual. These basis functions are system independent and easily computable, but the results 
are not always convincing and the number of necessary basis functions is almost unbeat- 
able. (Theoretically, plane waves could also be used for describing molecules, since the 
two-electron integrals and the expectation values are connected to the Fourier transform, 
thus they are easily computable, and this could balance the large number of necessary basis 
functions.) The reason, why so many plane waves are needed is that the wave functions 
around the nuclei need very high frequency terms, i.e. high resolution level, for repro- 
ducing the nuclear cusps. In the framework of Fourier analysis, the whole space has to 
be expanded at the same resolution, despite that in most of the space low frequency terms 
would be sufficient. 

Fully-numerical "basis-set free" Hartree-Fock (HF) calculations of atoms have been 
known since the 1960s (Vol. 1, pp. 322-326 and Vol. 2, pp. 15-30 of Ref. 03 and 
Refs. BUI HTl l42l 1431 ) and have proven helpful in constructing efficient finite basis sets 
for molecular calculations. In the late 1980s, Axel Becke used a fully-numerical density- 
functional theory (DFT) program for diatomics to show that many of the problems of DFT 
calculations at that time were due not to the functionals used, but rather numerical artifacts 
of the DFT programs of the 1970s EH .) Since that time, fully-numerical DFT codes have 
been implemented for polyatomic molecules using the finite element method (FEM), with 
PARSEC from the chemists point of view or OCTOPUS from the view of physicsts being 
a notable example. 

BlGDFT the pseudo potential code for bigger systems based as it is on traditional 
Hohenberg-Kohn-Sham DFT |[26l l27ll . could only calculate ground-state properties with 
an eye to order-N DFT. As a step to increase the feasibility of the code we formulated 
the wavelet-based linear-response time-dependent density-functional theory (TD-DFT) and 
here we support our first implementation for calculating electronic excitation spectra B31 . 
Electronic excitation spectra can be calculated from TD-DFT l46l using time-dependent 
linear response (LR) theory B71 l48l . Casida formulated LR-TD-DFT (often just refered 
to as TD-DFT) so as to resemble the linear-response time-dependent HF equations already 
familiar to quantum chemists RBI . That method was then rapidly implemented in a large 
number of electronic structure codes in quantum chemistry, beginning with the DEMON 
family of programs H91 and the TurboMol program ||50"1 . Among the programs that im- 
plemented "Casida's equations" early on was the FEM DFT program PARSEC lIBTI and 
also be found in the FEM DFT program OCTOPUS [52]]. See Ref. |53] for a recent FEM 
implementation of TD-DFT. Since a wavelet-based program offers certain advantages over 
these other FEM DFT programs, it was deemed important to also implement LR-TD-DFT 
in BlGDFT. 

In the next section we give a detailed description of the idea behind the multiresolution 
analysis and wavelets, with a historical note. Sec. [3] and Sec. [43 briefly presenting the the- 
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oretical introduction to DFT and TD-DFT, and Sec. [53 talks about the well-known Krylov 
space methods for solving eigenvalue equations involved in our implementation. Sec. [6] 
and Sec. [0 gives the numerical implementation of DFT and how we have implemented 
TD-DFT from the aspects of theoretical and algorithmic point of view on wavelets based 
pseudopotential electronic structure code BlGDFT, and in Section [8] we give the results 
of detailed comparisons between TD-DFT excitation spectra calculated with BlGDFT and 
with the implementation of Casida's equations in the GTO-based program deMon2k. The 
conclusion were drawn for future applications in the field of chemistry and some of the 
other problems are reviewed to draw chemists' greater attention to wavelets and to gain 
more benifits from using wavelet technique. 

2. Wavelet Theory 

The mathematics of wavelets is a fairly new technique, it can generally be used where 
one traditionally uses Fourier techniques. They incorporate the feature of having multiple 
scales, so very different resolutions can be used in different parts of space in a mathemat- 
ically rigorous manner. This matches many systems in nature well, for example molecule 
where the atomic orbitals are very detailed close to the cores, while they only vary slowly 
between them. Wavelet analysis can quite generally be viewed as a local Fourier analysis. 
From the wavelet expansion, or wavelet spectrum, of a function, /, it can be inferred not 
only how fast / varies, i.e. which frequencies it contains, but also where in space a given 
frequency is located. This property has important applications in both data compression, 
signal/image processing and noise reduction ll54l . Wavelet methods are also employed for 
solving partial differential equations |[55ll56l . and in relation to electronic structure methods 
a complete DFT program based on interpolating wavelets has been developed 11571 . 

2.1. The story of wavelets 

Most historical versions of wavelet theory however, despite their source's perspective, begin 
with Joseph Fourier. In 1807, a French mathematician, Joseph Fourier, discovered that all 
periodic functions could be expressed as a weighted sum of basic trigonometric functions. 
His ideas faced much criticism from Lagrange, Legendre and Laplace for lack of mathe- 
matical rigor and generality, and his papers were denied publication. It took Fourier over 
15 years to convince them and publish his results. Over the next 150 years his ideas were 
expanded and generalized for non-periodic functions and discrete time sequences. The fast 
Fourier transform algorithm, devised by Cooley and Tukey in 1965 placed the crown on 
Fourier transform, making it the king of all transforms. Since then Fourier transforms have 
been the most widely used, and often misused, mathematical tool in not only electrical en- 
gineering, but in many disciplines requiring function analysis. This crown however, was 
about to change hands. Following a remarkably similar history of development, the wavelet 
transform is rapidly gaining popularity and recognition. 

The first mention of wavelets was in a 1909 dissertation by Hungarian mathematician 
Alfred Haar. Haar's work was not necessarily about wavelets, as "wavelets" would not 
appear in their current form until the late 1980s. Specifically, Haar focused on orthogonal 
function systems, and proposed an orthogonal basis, now known as the Haar wavelet basis, 
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in which functions were to be transformed by two basis functions. One basis function is 
constant on a fixed interval, and is known as the scaling function. The other basis function 
is a step function that contains exactly one zero-crossing (vanishing moment) over a fixed 
interval (more on this later). 

The next major contribution to wavelet theory was from a 1930s French scientist Paul 
Pierre Levy. More correctly, Levy's contribution was less of a contribution and more of 
a validation. While studying the ins and outs of Brownian motion in the years following 
Haar's publication, Levy discovered that a scale-varied Haar basis produced a more accu- 
rate representation of Brownian motion than did the Fourier basis. Levy, being more of a 
physicist than mathematician, moved on to make large contributions to our understanding 
of stochastic processes. 

Contributions to wavelet theory between the 1930s and 1970s were slight. Most impor- 
tantly, the windowed Fourier transform was developed, with the largest contribution being 
made by another Hungarian named Dennis Gabor. The next major advancement in wavelet 
theory is considered to be that of Jean Morlet in the late 1970s. 

Morlet, a French geophysicist working with windowed Fourier transforms, discovered 
that fixing frequency and stretching or compressing (scaling) the time window was a more 
useful approach than varying frequency and fixing scale. Furthermore, these windows were 
all generated by dilation or compression of a prototype Gaussian. These window functions 
had compact support both in time and in frequency (since the Fourier transform of a Gaus- 
sian is also a Gaussian.) Due to the small and oscillatory nature of these window functions, 
Morlet named his functions as "wavelets of constant shape". In 1981, Morlet worked with 
Croatian-French physicist Alex Grossman on the idea that a function could be transformed 
by a wavelet basis and transformed back without loss of information, thereby outlining the 
wavelet transformation. It is of note that Morlet initially developed his ideas with nothing 
more than a handheld calculator. 

In 1986, Stephane Mallat noticed a publication by Yves Meyer that built on the concepts 
of Morlet and Grossman. Mallat sought Meyer's consult, and the result of said consult was 
Mallat's publication of multiresolution analysis. Mallat's MRA connected wavelet trans- 
formations with the field of digital signal processing. Specifically, Mallat developed the 
wavelet transformation as a multiresolution approximation produced by a pair of digital fil- 
ters. The scaling and wavelet functions that constitute a wavelet basis are represented by 
a pair of finite impulse response filters, and the wavelet transformation is computed as the 
convolution of these filters with the input function. The importance of Mallat's contribu- 
tion cannot be overstated. Without the fast computational means of wavelet transformation 
provided by the MRA, wavelets, undoubtedly, would not be the effective and widely used 
signal processing tools that they are today. 

In 1988, a student of Alex Grossman, named Ingrid Daubechies, combined the ideas of 
Morlet, Grossman, Mallat, and Meyer by developing the first family of wavelets as they are 
known today. Named the Daubechies wavelets, the family consists of 8 separate wavelet and 
scaling functions (more on this later). With the development of pair Daubechies wavelet and 
scaling functions is orthogonal, continuous, regular, and compactly supported, the founda- 
tions of the modern wavelet theory were laid. The last ten years mostly witnessed a search 
for other wavelets with different properties and modifications of the MRA algorithms. In 
1992, Albert Cohen, Jean Feauveau and Daubechies constructed the compactly supported 
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biorthogonal wavelets, which are preferred by many researchers over the orthonormal basis 
functions, whereas R. Coifman, Meyer and Victor Wickerhauser developed wavelet pack- 
ers, a natural extension of MRA. 

2.2. Multiresolution analysis 

A suitable gateway to the theory of wavelets is through the idea of MRA. A detailed de- 
scription of MRAs can be found in Keinert [58], from which a brief summary of the key 
issues are given in the following. 

A multiresolution analysis is an infinite nested sequence of subspaces L 2 (M) 

V® C Vj C ... C vp C ... (1) 

with the following properties 

• V°° is dense in L 2 

• f(x) G VJ 1 f(2x) G V™ +1 < n < oo 

• f{x) G Vp ^ f{x - 2~ n l) G VpO < I < (2 n - 1) 

• There exists a function vector ip of length j ; + 1 in L 2 such that 

{(Pj(x):0<k<j} 

forms a basis for V® . 

This means that if we can construct a basis of V®, which consists of only j + 1 functions, 
we can construct a basis of any space Vp, by simple compression (by a factor of 2 n ), and 
translations (to all grid points at scale n), of the original j + 1 functions, and by increasing 
the scale n, we are approaching a complete basis of L 2 . Since Vp C Vp +l the basis 
functions of Vp can be expanded in the basis of Vp +l 

tfix) d f 2"/V(2"x - = E fc ( W) ■ (2) 

i 

where h^s are the so-called filter matrix that describes the transformation between different 
spaces Vp. 

The MRA is called orthogonal if 

(<p$(x),<fi(x))=8 0l I j+1 , (3) 

where is the (j + 1) x (j + 1) unit matrix, and j + 1 is the length of the function 
vector. The orthogonality condition means that the functions are orthogonal both within 
one function vector and through all possible translations on one scale, but not through the 
different scales. 

Complementary to the nested sequence of subspaces Vp, we can define another series 
of spaces W™ that complements Vp in Vp +1 

yn+l = yn w n (4) 
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where there exists another function vector <p of length j + 1 that, with all its translations on 
scale n forms a basis for W^. Analogously to Eq. (0 the function vector can be expanded 
in the basis of V^ +1 

ff(x) d f 2 n / 2 <P(2 n x - I) = Y^9 {l) W + \x) . (5) 

l 

with filter matrices g( l \ The functions (f> also fulfill the same orthogonality condition as 
Eq. ©, and if we combine Eq. £T|) and Eq. we see that they must be orthogonal with 
respect to different scales. Using Eq. (0]) recursively we obtain 

V? = V} ®W$®W}®...® W^ 1 . (6) 
which will prove to be an important relation. 

2.3. Wavelets 

There are many ways to choose the basis functions ip and <fi (which define the spanned 
spaces VJ 1 and WJ 1 ), and there have been constructed functions with a variety of properties, 
and we should choose the wavelet family that best suits the needs of the problem we are 
trying to solve. (Wavelets are often denoted by %jj in the literature but the choice has been 
made here to denote them by so as to avoid confusion with the Kohn-Sham orbitals.) 
Otherwise, we could start from scratch and construct the new family, one that is custum- 
made for the problem at hand. Of course, this is not a trivial task, and it might prove more 
efficient to use an existing family, even though its properties are not right on cue. 

There is a one-to-one correspondence between the basis functions ip and cj>, and the filter 
matrices hS l > and gW used in the two-scale relation equations Eq. Q and Eq. ©, and most 
well-known wavelet families are defined only through their filter coefficients. 

In the following we are taking a different, more intuitive approach, for defining the 
scaling space V™ as the space of piecewise polynomial functions 

def 

V"' _ {/ : all polynomials of degree < j 
on the interval(2- rt Z, 2~ n {l + 1)) 

forO < I < 2 n , / vanishes elsewhere} . (7) 

It is quite obvious that one polynomial of degree j on the interval [0, 1] can be exactly 
reproduced by two polynomials of degree j, one on the interval [0, \] and the other on 
the interval [|, 1]. The spaces VJ 1 hence fulfills the MR A condition Eq. fl}, and if the 
polynomial basis is chosen to be orthogonal, the VJ 1 constitutes an orthogonal basis. 

2.4. An example: Simple Haar wavelets 

The basic wavelet ideas that we need can be easily explained using Haar wavelets |59l . 
These are simply the box functions shown in Fig. |2] We begin with a compact "mother 
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Figure 1. Wavelets (bottom) and scaling function (top). 

scaling function," in this case the Haar function, 

{0 ; x > 1 
1 ; < x < 1 . (8) 
; x < 

Translations, {<fi{x) = ip(x — i)}, of this mother function produces a crude basis set. Its 
relation to the grid of integers is obvious. Successively more refined basis sets may be 
generated by repeated application of the scaling operation consisting of contracting the 
functions to half their size in the x direction. The kth generation of scaling function is 
given by |<£>^(x) = (p(2 k x — Each generation has a fixed resolution related to an 
underlying grid with the same resolution. Let us now try to construct a multiresolution basis 
set. This is accomplished by (say) beginning with the third generation wavelets and taking 
sums and differences of adjacent functions until the eight third generation scaling functions 
have been replaced with the eight wavelets shown in Fig. [3] Notice how each generation of 
daughter wavelets is related to the mother wavelet by scaling, (jy*' (x) = <fr(2 k x — i). Notice 
also how the mother and two generations of daughter wavelets plus the mother scaling 
function (occasionally refered to as the "father wavelet") constitute a multiresolution basis 
set equivalent to the original third generation scaling basis set. Thus an arbitrary function, 
f(x), expressible in the original scaling basis, 

8 

/(*) = J>S 3) (*K (3) , (9) 
i=l 

has the wavelet transform, 

f (x) = ^(z)4°> + $Hz)4*> + £ ^(*k (1) + £ tfH X )4*> . do) 

i=0,l i=0,3 

Since the basis set is multiresolution, we may choose to add more grid points in some region 
of space and go locally to higher order wavelet expansions. It is also not always necessary 
to carry out a full wavelet transform, but rather it may be useful to just carry out a partial 
transform giving a linear combination of wavelets with several scaling functions at a time. 
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Figure 2. Haar scaling functions. 
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The extension to three dimensions is accomplished by using products of one-dimensional 
scaling functions and wavelets. Haar wavelets are just one type of wavelet basis set. It 
happens to be pedagogically useful but is not particularly useful for computations. 

2.5. Wavelet Basis 

The wavelet space W™ is defined, according to Eq. (@]), as the orthogonal complement of VJ 1 
in Vp +1 . The wavelet basis functions of WJ 1 are hence piece- wise polynomials of degree 
< j on each of the two intervals on scale n + 1 that overlaps with one interval on scale n. 
These piece-wise polynomials are then made orthogonal to a basis of VJ 1 and to each other. 
The construction of the wavelet basis follows exactly |6Ql where a simple Gram-Schmidt 
orthogonalization were employed to construct a basis that met the necessary orthogonality 
conditions. 

One important property of the wavelet basis is the number of vanishing moments. The 
j-th continuous moment of a function 4> is defined as the integrals 

dcf f 

fij _ / x J cj)(x)dx, (11) 
_ Jo 

and the function <p has M vanishing moments if 

p j= 0, fc = 0,...,M-l 

The vanishing momenets of the wavelet functions gives information on the approxima- 
tion order of the scaling functions. If the wavelet function <ft has M vanishing moments, 
any polynomial of order < M — 1 can be exactly reproduced by the scaling function tp, 
and the error in representing an arbitrary function in the scaling basis is of M-th order. By 
construction, x % is in the space Vj for < i < j, and since W® _L Vj, the first k + 1 
moments of <ff- must vanish. 

2.6. The scaling basis 

The construction of the scaling functions is quite straightforward, j+1 suitable polynomials 
are chosen to span any polynomial of degree j on the unit interval. The total basis for 
VJ 1 is then obtained by appropriate dilation and translation of these functions. Of course, 
any polynomial basis can be used, the simplest of them the standard basis {l,x, ...,x^}. 
However, this basis is not orthogonal on the unit interval. In the following, two choices of 
orthogonal scaling functions will be presented, and even though they span exactly the same 
spaces VJ 1 there are some important numerical differences between the two. 

In order to construct a set of orthogonal polynomials we could proceed in the same man- 
ner as for the wavelet functions and do a Gram-Schmidt orthogonalization of the standard 
basis {1, x, ...a? 7 }. If this is done on the interval x G [—1, 1] we end up with the Legendre 
polynomials {Lk} J k=0 - These functions are usually normalized such that Lfc(l) = 1 for all 
j. To make the Legendre scaling functions ip% we transform the Legendre polynomials to 
the interval x G [0, 1], and I? normalize 



ip%(x) = V2k + lL k (2x - 1), x € [0, 1] 



(12) 
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The basis for the space VJ 1 is then made by proper dilation and translation of ip k . Alpert et 
al. |[60l presented an alternative set of scaling functions with interpolating properties. These 
interpolating scaling functions ip k are based on the Legendre scaling functions {f^y k=0 , 
and the roots {yk}k=o an d weights {wk} J k=0 of the Gauss-Legendre quadrature of order 
j + I, and are organized in the linear combinations 

jp 

ip{(x) = ^^2^(y k )^(x), a: €[0,1], (13) 
i=0 

Again the basis of VJ 1 is made by dilation and translation of (f> k . The construction of (p\ 
gives them the interpolating property 

04, 

which will prove important for numerical efficiency. 

A detailed discussion on the properties of interpolating wavelets can be found in 
Donoho lloTll . but the case of interpolating wavelets is somewhat different. An important 
property of interpolating wavelets is the smoothness of any function represented in this 
basis. This property stems from general Lagrange interpolation. In the wavelet case the 
interpolating property applies within one scaling function vector only, which means that 
functions represented in this basis can be discontinous in any merging point between the 
different translations on any scale. 



2.7. Interpolating scaling functions 

Since the general introduction to wavelets were already made, we will now concencentrate 
our description on the level 3 interpolating scaling function (ISF) introduced by Deslauriers 
and Dubuc, and described in detail in Ref. (62]. Its main advantage is that it is fast and 
easy to perform nonlinear operations on functions represented in this basis, as long as the 
operation is local in shape. It also represents 3rd order polynomials exactly which means 
that it behaves very smoothly. 

We introduced the projection operator P n that projects an arbitrary function f(x) onto 
the basis Wji} of the scaling space V n (in the remaining of this text the subscript k of 
the scaling and wavelet spaces will be omitted, and it will always be assumed that we are 
dealing with a kth order polynomial basis.) 

f(x) « P n f(x) d f /"(*) = £ 5>jAfr(*) , (15) 

1=0 j=0 

where the expansion coefficients s™'/ , the so-called scaling coefficients, are obtained by 
the usual integral 

drf </. = f m^x)** , (i6) 

— Jo 

If this approximation turns out to be too crude, we double our basis set by increasing the 
scale and perform the projection p n + l _ This can be continued until we reach a scale N 
where we are satisfied with the overall accuracy of f N relative to the true function /. 



Wavelets for DFT and Post-DFT Calculations 



15 



In a perfect world, the projection in Eq. (fT6l ) could be done exactly, and the accuracy of 
the projection would be independent of the choice of polynomial basis. In the real world the 
projections are done with Gauss-Legendre quadrature and the expansion coefficients s™f 
of f(x) are obtained as 

2 _n(I+l) 

/ /(x^l^dx 

J2- nl 

2 -n/2 f l f { 2- n (x + l))^ (x)dx 

Jo 

kg-1 

2- n ' 2 Y,w q f{2- n {y q + l))$ fi {y q ) (17) 

q=0 

where {w q } q 'L are the weights and {yg} g l tne roots °f the Legendre polynomial L^ q 
used in k q ih order quadrature. 

By approximating this integral by quadrature we will of course not obtain the exact 
expansion coefficients. However, it would be nice if we could obtain the exact coefficients 
whenever our basis is flexible enough to reproduce the function exactly, that is if f(x) is 
a polynomial of degree < k. The Legendre quadrature holds a (2k — 1) rule which states 
that the fc-order quadrature is exact whenever the integrand is a polynomial of order 2k — 1. 
By choosing k q = k + 1 order quadrature we will obtain the exact coefficient whenever 
f(x) is a polynomail of degree < (k + 1) when projecting on the basis of fc-order Legendre 
polynomials. 

In the multidimensional case the expansion coefficients are given by multidimensional 
quadrature 

nj = 2 - dn/2 £ £ f(2- n (y q + l))Uf =1 w qil pl (y qi ) , (18) 

gi=0g2=0 q d =0 

using the following notation for the vector of quadrature roots 

Vq ^ (y<zi>2/? 2 >•••>%*) ) ( 19 ) 




This quadrature is not very efficient in multiple dimensions since the number of terms scales 
as (k + l) d . However, if the function / is separable and can be written f(x\,X2,---,x c i) = 
fi{xi)h{x2)--fd(xd), Eq. CEU can be simplified to 

k 

nf = 2 -dn/2 ud=i f t (2- n (y qt + k))w qi <p% >0 (y qi ) , (20) 

Qi=0 

which is a product of small summations and scales only as d(k + 1). 

The Legendre polynomials show very good convergence for polynomial functions f(x), 
and are likely to give more accurate projections. However, most interesting functions f(x) 
are not simple polynomials, and the accuracy of the Legendre scaling functions versus a 
general polynomial basis might not be very different. 
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By choosing the quadrature order to be k+ 1 a very important property of the Interpolat- 
ing scaling functions emerges, stemming from the specific construction of these functions 
Eq. (fT3T ). and the use of the k + 1 order quadrature roots and weights. The interpolating 
property Eq. (fPfl) inserts a Kronecker delta whenever the scaling function is evaluated in a 
quadrature root, which is exactly the case in the quadrature sum. This reduces Eq. (fTTT ) to 

sf = ^f(2- n ( Xj + I)) , (21) 
v i 

which obviously makes the projection k + 1 times more efficient. 

In multiple dimensions this property becomes even more important, since it effectively 
removes all the nested summations in Eq. ( fT8l ) and leaves only one term in the projection 

o-n/2 

S ;/ = /(2-' i (y J +0)nf =1 ^, (22) 

This means that in the Interpolating basis the projection is equally effective regardless of 
the separability of the function /. 



3. Density Functional Theory 

A method to resolve the electronic structure is by using variational principle 

Where {^>\H\^>) = f dr^*(r)H^(r), ^ denotes the electronic wavefunction and H the 
Hamiltonian. The energy computed from a guess ^ is an upper bound to the true ground 
state energy Eq. Full minimization of the functional E[*f>] will give the true ground state 
^f gs and energy E = E[^ 9S ]. 

Density-functional theory states that the many electron problem can be replaced by an 
equivalent set of self-consistent one-electron equations, the Kohn-Sham equations 

hf{r) = (~V 2 + v pp (r) + v H (r) + fi£ c (r)) ^f(r) = ef<(r) . (24) 

The eigenfunctions ipf are the one-electron wavefunctions that correspond to the minimum 
of the Kohn-Sham energy functional. In these wavefunctions, i is the orbital index and a 
denotes the spin, which can be either up f or down \. (spin a or (3.) 

The Hamiltonian H consists of four different parts: a part related to the kinetic energy 
of the electrons, the pseudopotentials v psp , the Hartree potential vjj and the exchange cor- 
relation potential v xc . The interaction of the positively charged nuclei with the electrons is 
described using the pseudopotential v psp instead of using the full Coulombic potential. The 
pseudopotential usually consists of both a local and a non-local part 

v psp (r) = v loc (r) + \l)*l(r, r ')( l \ ■ (25) 
l 
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The Hartree potential vh describes the interaction between electrons and is given by 

v H (r) = / dr -!— : f . (26) 



Finally, the exchange correlation potential v xc describes the nonclassical interaction be- 
tween the electrons and is given by the functional derivative of an exchange correlation 
energy functional 

= 6E f P J\ Pl) ■ ^ 
In these equations p a is the electron spin density, defined as 

p CT (r)=]T<|C(r)| 2 , (28) 

i 

where n° is the occupation number, i.e. the number of electrons in orbital i. In case of LDA 
(which we use throughtout this chapter) where there is no longer a distinction between spin 
up and spin down, orbitals can contain at most two electrons. 



4. Time-Dependent Density Functional Theory 

This section contains a brief review of the basic formalism of TD-DFT which is already 
well-known from the literature Roll . The time-dependent single particle Kohn-Sham equa- 
tions are, 

^_I V 2 + <; e// [p](r,i)J ^(r,*) =»J^Mr,t) (29) 
Here, the wave functions ijji(r,t) and v e ff[p](r,t) explicitly depend on time, whereas, 

v eff [p](r,t) = J2 v ion(r-'Ra)+ f y^lrdr' + v xc [p](r,t). (30) 

a J \ \ 

Using the adiabatic approximation (AA), (which is local in time) 

5E xc [p] 



v xc [p](r,t) 



Sp(r) 



5v xc [p](r,t) _ S 2 E xc [p] 
dp(r',t) ~ 1 j ,5p(r)Mr')' { ' 

and using the LDA, 

E xc [p] = J p(r)e xc (p(r))dr . (32) 



The method that we wish to use here is Casida's approach ||48l . This section explains 
the deriving equations of linear-response (LR) TD-LDA method. 

The time-dependent perturbation to the external potential can be written as, 



$ v eff[p]{r,t) = 5v app i(r,t) + 6v S cF[p](r,t) 



(33) 
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where, 



vsCF[p](r,t) 



— dr + VxcHlr ,t) 



r — r 



(34) 



The LR of the density matrix (DM) can be written in terms of generalised susceptibility \ 
as below, 

5P ija {uj) =^2xija,klr(^ v tU( U} ) ( 35 ) 
klr 

After some mathematical steps, one can end-up with the sum-over-states (SOS) representa- 
tion of x, 



(36) 



W - (uJlkr) 

where A^ T = n^ T — n kr the difference in occupation numbers and w^. T = e kr — e; T 
the difference between the eigenvalues of Ith and kth states. In the basis of Kohn-Sham 
orbitalsj^jcr}, we can re-write the LR-DM equation as, 



5P ija (u) 



A 



(37) 



Now the term 5v SCF is complicated because it itself depends on the response of the DM. 

5vf^ F (uj) = K ijaMr 5P kl r(oj) (38) 

Where, 



K 



ija,klr 



dv SCF 



dP 



klr 



whose integral form is, 



K 



ijcr,klT 



V40)Vv( r ) 



r — r 



6 2 E xc [p] 
Sp a (r)Sp T (r') 



(39) 



i; kT (r')tfr(r')drdr' (40) 



If the response is due to a real spin independent external perturbation, Sv appl , then we can 
restrict ourselves to the real density response and the coupling matrix is symmetric. 
After some algebra, the real parts of the DM elements R5P(lo) can be given as, 



klr 



^klr^klr 



(id — oj Ut ) — 2Kjj (7)HT 



R(5P klT )(u) = 5v^ l (u) (41) 



Here the real part of RSP a (uj) means the Fourier transform of the real part of $t5P a (t). 
Thus the real part of the first-order DM obtained from the solution of the above linear equa- 
tions gives access to the frequency-dependent polarizabilities. This leads to the following 
eigenvalue equation from which the excitation energies and oscillator strengths can be ob- 
tained. 



QFj = ojjFj , 



where, 



rtij<7,klr — fiikdjlficrT^klT + ^ij(T UJ ija'^-ijcr,klT\/ ^klr^klr 



(42) 
(43) 
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Here, the desired excitation energies are equal to oj\ and the oscillator strengths // are ob- 
tained from the eigenvectors Fj . The frequency-dependent polarizability is directly real ted 
to vertical excitation energies, oscillator strength and transition dipole moments /i/, 



^— ' loj z — ur 3 coi z — uj z 



l 



5. Krylov Space Methods 

The methods described in this article involve solving very large eigenvalue problems. One 
of these is the matrix form of the Kohn-Sham orbital equation Eq. (l24l while the other is the 
LR-TD-DFT equation Eq. (l42l ). The first is very large because the wavelet basis set is very 
large while the other is very large because it is of the order of the number of unoccupied 
orbitals times the number of unoccupied orbitals on each side. It is important to realize that 
special methods must be and are being used to solve these very large eigenvalue problems. 
In particular, BlGDFT make use of the block Davidson variant of the Krylov space method 
to solve the Kohn-Sham equation while BlGDFT and most other codes solving the LR-TD- 
DFT equation Eq. (|42l) also make use of the the block Davidson method. Krylov methods 
and the block Davidson method are briefly described in this section. 

Krylov space methods may be traced back to a paper in the early 1930s written by 
the Russian mathematician Alexei Nikolaevich Krylov. The main idea is that to solve a 
matrix problem involving a matrix A, it is frequently never actually necessary to construct 
the matrix A because iterative solutions only require a reasonable first guess followed by 
repeated action of the operator iona vector. A number of such methods are known with 
Lanczos diagonalization and the discrete inversion in the iterative subspace (DIIS) ll63l as 
particlarly well-known examples in theoretical chemical physics. Given a vector x, the 
Krylov space of order r is given by, 

/C r (A, x) = span {x, Ax, A 2 x, ■■■ , A r x} . (45) 

The Davidson diagonalization method ll64l for solving the matrix eigenvalue problem 

Ax = ax , (46) 



is deceptively simple. Suppose that we want the lowest eigenvalue and eigenvector and we 
have an intial guess vector, x^ . Then we can always write, 

x = x w +5x, (47) 



is the component of the exact solution which is orthogonal to the intial gues vector. Sim- 
ple algebra then gives a formula highly reminiscent of Rayleigh-Schrodinger perturbation 
theory but exact, 

5x = [Q(al- A)Q]~ 1 (A-al)f (0) , (48) 

where, 

Q = 1 - f(% (°)t , (49) 
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projects onto the subspace orthogonal to the guess vector. Solving Eq. (|48T ) requires us to 
overcome two difficulties. The first is that we need a guess for the eigenvalue a, but this 
is easily remedied by taking = f' ^Ax' ' /x^x^ and then iterating. The problem 
greater problem is to invert the matrix [Q(al — A)Q]. It actually turns out that a highly 
accurate inversion is not really needed (and actually can cause some problems.) Instead it 
is better just to replace this difficult inversion with, 

5x = (al-D)~ 1 (A-al)f {0) , (50) 

where D is some diagonal matrix, hence easy to invert. Orthogonalizing 5x to x^ ' and 
normalizing produces x^\ which is the next basis vector in our Krylov space. Setting up 
and diagonalizing the 2x2 matrix of A in this basis and taking the lowest eigenvalue gives 
us the next estimate a^ l \ If application of Eq. (|48T ) is close to zero then we have solved 
the eigenvalue problem Eq. (l46l) . otherwise we have a new Sx with which to generate x^ 
and so on and so forth until convergence. The block Davidson method ll65l extends the 
Davidson method to the lowest several eigenvalues and eigenvectors. 

Davidson diagonalization works well when started from a reasonably good initial guess, 
otherwise the Lanczos method may be advantageous. One of the most recent incarnations of 
the Lanczos method is the Liouville-Lanczos method for solving the LR-TD-DFT problem 
MMMM- 

6. Numerical Implementation of DFT in BigDFT 

Computational physics/chemistry is the transformation and implementation of scientific 
theory into efficient algorithms which requires both theoretical and experimental skill. The 
transformation of a new theory into an efficient algorithm requires understanding of pro- 
gramming concepts, mathematical and physical intuition and theoretical insight, whereas 
the production of the computer code is much like experimentation, requiring debugging, 
testing and organisation to yield a highly efficient product. It is also an adaptation of new 
scientific theory into computer code exploiting the advances in compiler, programming lan- 
guage and hardware technology. The aim is to afford an algorithm to enable efficient com- 
putation, portability of code, ease of adaptability and to document the science. To afford 
such an algorithm requires an intuitive understanding of the physics to be implemented, 
much experimentation with optimisation and debugging of the developing code, a suitable 
choice of programming language, as well as a basic overview of the nature of the platforms 
for which the code is intended. 

The Kohn-Sham scheme of DFT greatly reduces the complexity of ground state elec- 
tronic structure calculations by recasting the many-body problem into a (self-consistent) 
single-particle problem. For real atomistic systems, however, the KS equations are still 
difficult to solve and further approximate techniques are required. In general it is impor- 
tant, though, that these approximations can be controlled in such a way that the associated 
error does not exceed the error introduced by the xc-functional. While DFT accounts for 
approximately 90% of all quantum chemical calculations being performed, the sometimes 
unpredictable nature of results and the inability to systematically improve the quality of 
calculation may mean that a place for the conventional correlated techniques remains in the 
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quantum chemist's tool kit. In this work the detailed description of DFT program BlGDFT 
has been given. BlGDFT 11571 has been developed as an European project (FP6-NEST) 
from 2005 to 2008, and is a wavelet-based pseudopotential implementation of DFT and 
TD-DFT. For complimentary purpose, Gaussian based quantum chemistry DFT code DE- 
Mon2k f70l is also used but we are not going to discuss the numerical implementation of 
deMon2k here and we restrict ourselves to use deMon2k for revalidating our recent im- 
plementation of LR-TD-DFT in BlGDFT. However in the following sections, we are going 
to recast how the fundamental computational operations were performed in BlGDFT. 

6.1. Daubechies Wavelets 

Before embarking on our own endeavours, we should make some reference to related work. 
First, it should be acknowledged that a considerable amount of work has been done already 
in pursuit of a wavelets in the electronic structure calculations fTTl 1721 1731 1741 1751 1761 . 
The object of using wavelets as basis set is to associate an expansion coefficient to each 
of the piece-wise wavelets. The expansion coefficients are free to vary from one wavelet 
function to the next. This feature enables wavelets as highly localised continuous functions 
of a fractal nature that have finite supports. The Daubechies wavelets have no available 
analytic forms, and they are not readily available in sampled versions. They are defined 
effectively by the associated dilation coefficients. These express a wavelet in high resolution 
and a scaling function in the low resolution-which has the same width and which stretches 
to zero-as a linear combination of the more densely packed and less dispersed scaling 
functions that form a basis for the two resolution level in combination. 

The fact that the Daubechies wavelets are known only via their dilation coefficients is no 
implediment to the discrete wavelet transform. This transform generates the expansion co- 
efficients associated with the wavelet decomposition of a data sequence. In this perspective, 
the dilation coefficients of the wavelets and of the associated scaling functions are nothing 
but the coefficients of a pair of quadrature mirror filters that are applied succesively. 

As described above Daubechies family consists of two fundamental functions: the scal- 
ing function 4>(x) and the wavelet ip(x) (see Fig. @]) The full basis set can be obtained from 
all translations by a certain grid spacing h of the scaling and wavelet functions centered at 
the origin. These functions satisfy the fundamental defining (refinement) equations, 

m 

4>(x) = V2 h^(2x-j), 

j=X—m 
in 

<p(x) = V2 9j<l>(2x-j). (51) 

j=l— m 

which relate the basis functions on a grid with spacing h and another one with spacing h/2. 
The coefficients, hj and gj , consitute the so-called "filters" which define the wavelet family 
of order m. These coefficients satisfy the relations, ^ ■ hj = 1 and gj = (— l) J 7i_j + i. 
Eq. (I5TI ) is very important since it means that a scaling-function basis defined over a fine 
grid of spacing h/2 may be replaced by combining a scaling-function basis over a coarse 
grid of spacing h with a wavelet basis defined over the fine grid of spacing h/2. This then 
gives us the liberty to begin with a coarse description in terms of scaling functions and 
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then add wavelets only where a more refined description is needed. In principle the refined 
wavelet description may be further refined by adding higher-order wavelets where needed. 
However in BlGDFT we restricted ourselves to just two levels: coarse and fine associated 
respectively with scaling functions and wavelets. 

For a three-dimensional description, the simplest basis set is obtained by a tensor prod- 
uct of one-dimensional basis functions. For a two resolution level description, the coarse 
degrees of freedom are expanded by a single three dimensional function, ^ ia 4 (r), while 
the fine degrees of freedom can be expressed by adding another seven basis functions, 
fiji 32 33 ( r )' wrucn include tensor products with one-dimensional wavelet functions. Thus, 
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Figure 4. Daubechies scaling function cj)(x) and wavelet ip(x) of order 16. 
the Kohn-Sham wave function ^(r) is of the form 

7 

hM,ia 31,32,33 "i 

The sum over h,i2,h runs over all the grid points contained in the low-resolution regions 
and the sum over J1J2J3 runs over all the points contained in the (generally smaller) 
high resolution regions. Each wave function is then described by a set of coefficients 
- 2 - 3 }, v = 0, ...,7. Only the nonzero scaling function and wavelet coefficients are 
stored. The data is thus compressed. The basis set being orthogonal, several operations 
such as scalar products among different orbitals and between orbitals and the projectors of 
the nonlocal pseudopotential can be directly carried out in this compressed form. 

6.2. Treatment of kinetic energy 

The matrix elements of the kinetic energy operator among the basis functions of our mixed 
representation (i.e., scaling functions with scaling functions, scaling function with wavelets 
and wavelets with wavelets) can be calculated analytically f77l . For simplicity, let us il- 
lustrate the application of the kinetic energy operator onto a wavefunction ip that is only 
expressed in terms of scaling functions. 

V>0; y,z) = Y s h,i2,i-A( x / h - h) 4>(y/h - i 2 ) 4>(z/h - i 3 ) 



Wavelets for DFT and Post-DFT Calculations 



23 



The result of the application of the kinetic energy operator on this wavefunction, projected 
to the original scaling function space, has the expansion coefficients 

*ii,i2,*3 = ~ 7^3 J <f>(x/h - h) 4>{y/h - i 2 ) 4>{z/h - i 3 ) x 

x V 2 i(j(x, y, z)dxdydz . 
Analytically the coefficients 8i u i a j 3 and Si lt i 2i i 3 are related by a convolution 

where 

Ki\,i2,iz = Ti 1 Ti 2 Ti 3 , (54) 
where the coefficients Tj can be calculated analytically via an eigenvalue equation: 

f d 2 

Ti = cj)(x)-Q-2<p(x - i)dx 

<j>(2x - v)—$(2x -2i- n)dx 
f d 2 

= ^ h u h^2 2 / 0(y)^-70(y -2i-[i + v)dy 

= ^2 huh^ 2 2 T2i- v + lL 

Using the refinement equation Eq. (I5TI ). the values of the Tj can be calculated analytically, 
from a suitable eigenvector of a matrix derived from the wavelet filters [77]] • For this reason 
the expression of the kinetic energy operator is exact in a given Daubechies basis. 

Since the 3-dimensional kinetic energy filter K^^^ is a product of three one- 
dimensional filters Eq. (l54l the convolution in Eq. (l53l can be evaluated with 3N1N2N3L 
operations for a three-dimensional grid of N1N2N3 grid points. L is the length of the 
one-dimensional filter which is 29 for our Daubechies family. The kinetic energy can thus 
be evaluated with linear scaling with respect to the number of nonvanishing expansion co- 
efficients of the wavefunction. This statement remains true for a mixed scaling function- 
wavelet basis where we have both nonvanishing s and d coefficients and for the case where 
the low and high resolution regions cover only parts of the cube of NiN 2 N s grid points. 



6.3. Treatment of local potential energy 

In spite of the striking advantages of Daubechies wavelets the initial exploration of this 
basis set f78l did not lead to any algorithm that would be useful for practical electronic 
structure calculations. This was due to the fact that an accurate evaluation of the local 
potential energy is difficult in a Daubechies wavelet basis. 
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By definition, the local potential v(r) can be easily known on the nodes of the uniform 
grid of the simulation box. Approximating a potential energy matrix element «ij fc;i' ,j>,k' 



Vij,k;i',j>,k' = J dr0 i / ij / ik /(r)v(r)(/) i j ik (r) 

by 

Vi,j,k;i',j',k' ~ ^2 4>i' ,j' ,k'{ r l,m,n) v { Y l,m,n)4>ij,k{ v l,m,n) 
l,m,n 

gives an extremely slow convergence rate with respect to the number of grid points used to 
approximate the integral because a single scaling function is not very smooth, i.e., it has a 
rather low number of continuous derivatives. A. Neelov and S. Goedecker IT791 have shown 
that one should not try to approximate a single matrix element as accurately as possible 
but that one should try instead to approximate directly the expectation value of the local 
potential. The reason for this strategy is that the wavefunction expressed in the Daubechy 
basis is smoother than a single Daubechies basis function. A single Daubechies scaling 
function of order 16 (i.e., the corresponding wavelet has 8 vanishing moments) has only 2 
continuous derivatives. More precisely its index of Holder continuity is about 2.7 and the 
Sobolev space regularity with respect to p = 2 is about 2.91 (80']. A single Daubechies 
scaling function of order 16 has only 4 continuous derivatives. By suitable linear com- 
binations of Daubechies 16 one can however exactly represent polynomials up to degree 
7, i.e., functions that have 7 non-vanishing continuous derivatives. The discontinuities get 
thus canceled by taking suitable linear combinations. Since we use pseudopotentials, our 
exact wavefunctions are analytic and can locally be represented by a Taylor series. We 
are thus approximating functions that are approximately polynomials of order 7 and the 
discontinuities nearly cancel. 

Instead of calculating the exact matrix elements we therefore use matrix elements with 
respect to a smoothed version <f> of the Daubechies scaling functions. 



Vi,j,k;i',j',k' ~ ^2 $V j\k'{ Y l,m,n) v { v l,m,n)$i,j,k{ Y l,m,n) 
l,m,n 

= ^0,0,o( r l-i\rri-j',n-k')y( r l,rn,n)(f>0,0,o( r l-i,rn-j,n-k) (55) 



l,m,n 

where the smoothed wave function is defined by 

<Po,0,o( r l,m,n) = UiLO m Ld n 



and bj\ is the "magic filter". Even though Eq. (1551 ) is not a particulary good approximation 
for a single matrix element it gives an excellent approximation for the expectation values of 
the local potential energy 



dx J dy J dz^(x,y,z)v(x,y,z)ip(x,y,z) 
and also for matrix elements between different wavefunctions 

dx dy dzil)i{x,y,z)v{x,y,z)il)j{x,y,z) 



Wavelets for DFT and Post-DFT Calculations 



25 



in case they are needed. Because of this remarkable achievement of the filter uj we call it 
the magic filter. 

Following the same guidelines as the kinetic energy filters, the smoothed real space 
values rpij t k of a wavefunction ifi are calculated by performing a product of three one- 
dimensional convolutions with the magic filters along the x, y and z directions. For the 
scaling function part of the wavefunction the corresponding formula is 



12,13 



E • 

h ,32,33 



S 3l,32,33 V i 1 



(1) (1) 

-2ji V i 2 - 



(1) 



-2j 3 



where v\ is the filter that maps a scaling function on a double resolution grid. Similar 
convolutions are needed for the wavelet part. The calculation is thus similar to the treatment 
of the Laplacian in the kinetic energy. 

Once we have calculated ipij t k the approximate expectation value ey of the local po- 
tential v for a wavefunction ip is obtained by simple summation on the double resolution 
real space grid: 

e v = ,32,33 V 3^,32,33^3^,32,33 
31,32,33 



6.4. Treatment of the non-local pseudopotential 

The energy contributions from the non-local pseudopotential have for each angular moment 
/ the form 

^2(ip\Pi)hij(pj\^} 

where \pi) is a pseudopotential projector. Once applying the Hamiltonian operator, the 
application of one projector on the wavefunctions requires the calculation of 

y 

If we use for the projectors the representation of Eq. (l52l (i.e., the same as for the wave- 
functions) both operations are trivial to perform. Because of the orthogonality of the basis 
set we just have to calculate scalar products among the coefficient vectors and to update the 
wavefunctions. The scaling function and wavelet expansion coefficients for the projectors 
are given by ||8T1 

ypw^i^^wd 1 " . y"p( r ) ( r ) dr • (56) 

The GTH-HGH pseudopotentials ll82l l83l have projectors which are written in terms 
of gaussians times polynomials. This form of projectors is particularly convenient to be 
expanded in the Daubechies basis. In other terms, since the general form of the projector is 

(r\p) = e~ cr2 x ix y iy z iz , 
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the 3-dimensional integrals can be calculated easily since they can be factorized into a 
product of 3 one-dimensional integrals. 



The one-dimensional integrals are calculated in the following way. We first calculate 
the scaling function expansion coefficients for scaling functions on a one-dimensional grid 
that is 16 times denser. The integration on this dense grid is done by the well-known quadra- 
ture introduced in ll84l . that coincides with the magic filter |f79l . This integration scheme 
based on the magic filter has a convergence rate of /i 16 and we gain therefore a factor of 
16 16 in accuracy by going to a denser grid. This means that the expansion coefficients are 
for reasonable grid spacings h accurate to machine precision. After having obtained the 
expansion coefficients with respect to the fine scaling functions we obtain the expansion 
coefficients with respect to the scaling functions and wavelets on the required resolution 
level by one-dimensional fast wavelet transformations. No accuracy is lost in the wavelet 
transforms and our representation of the projectors is therefore typically accurate to nearly 
machine precision. In order to treat with the same advantages other pseudopotentials which 
are not given under the form of gaussians it would be necessary to approximate them by a 
small number of gaussians. 

6.5. The Poisson operator 

Solving the Poisson equation for an arbitrary charge distribution is a non-trivial task, and 
is of major importance in many field of science, especially in the field of computational 
chemistry. A huge effort has been put into making efficient Poisson solvers, and usual 
real-space approaches includes finite difference (FD) and finite element (FE) methods. FD 
is a grid-based method, which is solving the equations iteratively on a discrete grid of 
pointvalues, while FE is expanding the solution in a basis set, usually by dividing space 
into cubic cells and allocate a polynomial basis to each cell. 

It is well-known fact that the electronic density in molecular systems is rapidly varying 
in the vicinity of the atomic nuclei, and a usual problem with real-space methods is that an 
accurate treatment of the system requires high resolution of gridpoints (FD) or cells (FE) in 
the nuclear regions. Keeping this high resolution uniformly throughout space would yield 
unnecessary high accuracy in the interatomic regions, and the solution of the Poisson equa- 
tion for molecular systems is demanding a multiresolution framework in order to achieve 
numerical efficiency. This chapter is concerned with a way of doing DFT and TD-DFT cal- 
culations, one where the multiresolution character is inherent in the theory, namely using 
wavelet bases. 

In order to evaluate the Hartree potential, we need to rewrite the standard Poisson equa- 
tion to an integral form. The equation, in its differential form, is given as 




(57) 



(58) 



V 2 -u(x) = 47T/?(x) , 



(59) 
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where p(x) is the known (charge) distribution, and u(x) is the unknown (electrostatic) 
potential. It is a standard textbook procedure to show that the solution can be written as the 
integral 

«(x) = J G(x,y)p(y)dy, (60) 

where G(x, y) is the Green's function which is the solution to the fundamental equation 
with homogeneous (Dirichlet) boundary conditions 

V 2 G(x,y) = 5(x-y) 

G(x, y) = 0,xG boundary (61) 

This equation can be solved analytically and the Green's function is given (in three dimen- 
sions) simply as 

G(x,y) = — ^ (62) 

ll x >y|| 

This is the well-known potential arising from a point charge located in the position y , which 
is exactly what Eq. (loTT) describes. 



6.6. Numerical separation of the kernel 

The Green's function kernel as it is given in Eq. (l62l is not separable in the cartesian co- 
ordinates. However, since we are working with finite precision we can get by with an 
approximate kernel as long as the error introduced with this approximation is less than our 
overall accuracy criterion. If we are able to obtain such a numerical separation of the kernel, 
the operator can be applied in one direction at the time, allowing us to use the expressions 
derived above for one-dimentional integral operators to solve the three-dimensional Pois- 
son equation. This is of great importance because it reduces the scaling behavior to become 
linear in the dimension of the system. 

The Poisson kernel can be made separable by expanding it as a sum of Gaussian func- 
tions, specifically 

1 M e 

-~^ Wfe e-^ 2 . (63) 

k=l 

where and p^ are parameters that needs to be determined, and the number of terms 
M € , called the separation rank, depends on the accuracy requirement and on what interval 
this expression needs to be valid. Details of how to obtain this expression can be found in 
ll7T1l72ll . and will not be treated here, but it should be mentioned that the separation rank is 
usually in the order of 100, e.g, it requires M e = 89 to reproduce £ on the interval [10 -9 , 1] 
in three dimensions with error less than e = 10~ 8 . 

Finally, figure [5] summarize this complete section into a flow-chart type diagram. This 
kind of explanation is necessary for begineers because there are different functions used for 
the different operations in BlGDFT. As one can see from the figure, The KS wavefunctions 

are expressed in terms of Daubechies wavelets and the projection of Hamiltonian V n i\ip) 
and of pseudopotential operators | Hip) also expressed using Daubechies wavelets. The rest 
of the operations such as kinetic energy \V 2 tp), potential energy operator V(x)ip(x), and 
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the local densities p(x) are all expressed using interpolating scaling functions, in which 
the Hartree Vh(x), local part potential energy Vi oc (x) and xc operations V xc (x) were per- 
formed in real space. The interconnecting lines between different operations represents the 
transformation between Daubechies wavelets-to-ISFs or the transformation of real space- 
to-fourier space representation. 



Daubecies 
wavelets 




Interpolating 
scaling functions 



Real space 




Figure 5. Operations performed in BlGDFT 



7. BigDFT and TD-DFT 



We want to solve Casida's equation [48], 

Y ah bh \ / 1 \] ( X 
\B*(u) A*(uj)J U {0 -l)\\Y 

where X and Y represents the pseudo eigenvectors; the matrices A and B are defined as 



0, 



(64) 



and, 



in which the integral form of the coupling matrix K is given by, 



K 



pqa,rsT 
1 



+ 



d 2 E xc [p] 



r — f| dp (7 {r)dp T {r l ) 



(65) 
(66) 



(67) 
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The universal adiabatic approximation is applied to Eq. (I67T ) to remove the frequency de- 
pendence of the kernel. 

The electronic transitions occur with an infinitesimal perturbation obtains the above 
described non-Hermitian eigenvalue Eq. (l64l) . Where the response is due to a real spin 
independent external perturbation, and the actual response is described as the real density 
response. However, an unitary transformation is necessary to convert Eq. (l64l into the real 
eiganvalue problem. In Eq. (l64l) . all occupied-occupied and virtual-virtual element contri- 
butions are zero whereas only the elements that are from virtual-occupied and occupied- 
virtual parts are taken into account. Moreover if we only restricted to virtual-occupied 
elements and neglecting the occupied-virtual elements of Eq. (l64l leads to a Hermitian 
eigenvalue equation of the dimension one-half of that TD-DFT working equation is said to 
be Tamm-Dancoff approximation (TDA) and it is written as, 

AX = wX, (68) 

where A is as same as in Eq. (|42l) . The matrix A is just restricted to number of single 
excitations. 

However, the explicit form of Eq. (l68l) is, 

n(w)J> = w 2 F/ , (69) 

where 

^iacrjbr — ^ia^jb^uri^aa €ia) (70) 
2 yi^acT £icr)Kiacr,jbT yi^aa ^ia) i 

where e; irj — e aa is the energy eigenvalue differences of i th and a th states. Solving Eqs. (l69l 
yields TD-DFT excitation energies to and F/'s are the corresponding oscillator strengths 
which are defined from the transition dipole moments. 



7.1. Calculation of Coupling Matrix 

We are now in a position to understand the construction of the coupling matrix Eq. (l67l) in 
our implementation of TD-DFT in BlGDFT, which we split into the Hartree andxc parts, 

K a iu,bjr = K a icr,bjT ^aja,bjr • (71) 

Instead of calculating the Hartree part of coupling matrix directly as, 

K* a , bjT = J jraa(r)Mr)y^^br(r'W n (r')drdr', (72) 
we express the coupling matrix element as, 

K aia,bjr = J <a Wi<r( r )%>( r ) , (73) 

where, 

vaia(v) = J r^4dr', (74) 
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and, 

Paia(r) = VC( r )?M r ) • ( 75 ) 

The advantage of doing this is that, although p a i a and v a i a are neither real physical charge 
densities nor real physical potentials, they still satisfy the Poisson equation, 

V 2 iw(r) = -47r/w(r) , (76) 

and we can make use of whichever of the efficient wavelet-based Poisson solvers already 
available in BlGDFT, is appropriate for the boundary conditions of our physical problem. 

Once the solution of Poisson's equation, v a i a (r), is known, we can then calculate the 
Hartree part of the kernel according to Eq. (l73l) . Inclusion of the xc kernel is accomplished 
by evaluating, 

Kaia,bjr = J M aia (r)p bjT (r) dr , (77) 

where, 

M ma {r) = v aia (r) + f p ai a(r')f^ T (r, r') dr' . (78) 

We note that fxi T {r, r') = f°c (r, r')5(r — r') for the LDA, so that no integral need actually 
be carried out in evaluating M a i a (r). The integral in Eq. (1771 ) is, of course, earned out 
numerically in practice as a discrete summation. 



8. Results 

We now wish to illustrate a bit how wavelet calculations work in the BlGDFT pro- 
gram. Comparison will be made against results obtained with the GTO-based program 
deMon2k. This work is very similar to our previous work reporting the first implemen- 
tation of wavelet-based TD-DFT with illustration for N2 and application to the absorption 
spectrum of a medium-sized organic molecule of potential biomedical use as a fluores- 
cent probe RSI . Here however we will present new BlGDFT results for a different small 
molecule, namely carbon monoxide. Though CO is roughly isoelectronic with N2, CO has 
the interesting feature of having a low-lying bright state in its absorption spectrum. 

8.1. Computational Details 

Calculations were carried out with deMon2k and BlGDFT with the LDA-optimized bond 
length of 1.129 A. 

8.1.1. deMon2k 

deMon2k resembles a typical GTO-based quantum chemistry program in that all the in- 
tegrals other than the xc-integrals, can be evaluated analytically. In particular, deMon2k 
has the important advantage that it accepts the popular GTO basis sets common in quan- 
tum chemistry and so can benefit from the experience in basis set construction of a large 
community built up over the past 50 years or so. In the following, we have chosen to use 
the well-known correlation-consistent basis sets for this study ll85l [86l . (Note, however, 
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that the correlation-consistent basis sets used in deMon2k lack f and g functions but are 
otherwise exactly the same as the usual ones.) The advantage of using these particular basis 
sets is that there is a clear hierarchy as to quality. 

An exception to the rule that integrals are evaluated analytically in deMon2k are the 
xc-integrals (for the xc-energy, xc -potential, and xc -kernel) which are evaluated numeri- 
cally over a Becke atom-centered grid. This is important because the relative simplicity of 
evaluating integrals over a grid has allowed the rapid implemenation of new functionals as 
they were introduced. We made use of the fine fixed grid in our calculations. 

As described so far, deMon2k should have 0(N 4 ) scaling because of the need to 
evaluate 4-center integrals. Instead deMon2k uses a second atom-centered auxiliary GTO 
basis to expand the charge density. This allows the the elimination of all 4-center integrals 
so that only 3-center integrals remain for a formal 0(N 3 ) scaling. In practice, integral 
prescreening leads to 0{NM) scaling where M is typically between 2 and 3. We made use 
of the A3 auxiliary basis set from the deMon2k automated auxiliary basis set library. 

All calculations were performed using standard deMon2k default criteria. The im- 
plementation of TD-DFT in deMon2k is described in Ref. [87]. (The charge density 
conservation constraint is no longer used in deMon2k TD-DFT calculations.) Although 
full TD-LDA calculations are possible with deMon2k, the TD-LDA calculations reported 
here all made use of the TDA. 



8.1.2. BlGDFT 
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(a) H2O in a box 



(b) H2O in a box showing fine grid (c) H2O in a box showing coarse 
resolution grid resolution 



Figure 6. Adaptive grid in BlGDFT |(a)[|(b)| and[(c) 



The main thing to vary in BlGDFT is the grid which is of more profound importance 
than in deMon2k because it is the grid which supports the wavelets. Figures [6(a)[ |6(b)[ 
and |6(c)| give an idea of what the grid looks like for the small familial - molecule of water. 
Conceptually the molecule is in a very large box (Fig. |6(a)| ) A fine grid is placed in the 
regions of high electron density around the molecule (Fig. |6(b)] ) A coarse grid is used in 
a larger region where the electron density varies more slowly (Fig. |6(c)[ ) The BlGDFT 
grid is characterized by the triple /i 9 /crmult/frmult. The first number in the triple (h g ) is 
a real number which specifies the nodes of the grid in atomic units. The second number 
(the integer-valued crmult) is the coarse grid multiplier. And the third number (the integer- 
valued frmult) is the fine grid multiplier. Two points must be clearly understood when 
looking at Figs. |6(a)||6(b)[ and |6(c)| The first is that, while the box may determine the limits 
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Table 1. Basis set dependence of the HOMO and LUMO energies and of the HOMO- 
LUMO gap (eV) calculated using deMon2k. 



Basis Set 


— (HOMO 


—clumo 


^HOMO-LUMO 


STO-3G 


-5.5350 


1.2428 


4.2922 


DZVP 


-8.9271 


-2.0942 


6.8329 


TZVP 


-9.0287 


-2.1902 


6.8385 


CC-PVDZ 


-8.6729 


-1.7823 


6.8906 


CC-PVTZ 


-9.0419 


-2.1195 


6.9224 


CC-PVQZ 


-9.0944 


-2.1971 


6.8973 


CC-PV5Z 


-9.1169 


-2.2400 


6.8769 


CC-PCVDZ 


-8.6905 


-1.7922 


6.8983 


CC-PCVQZ 


-9.0957 


-2.1988 


6.8969 


CC-PCVT7 


-9.0371 


-2.1165 


6.9206 


CC-PCV5Z 


-9.1172 


-2.2401 


6.8771 


AUG-CC-PVDZ 


-9.0910 


-2.2345 


6.8565 


AUG-CC-PVQZ 


-9.1286 


-2.2567 


6.8719 


AUG-CC-PVTZ 


-9.1306 


-2.2535 


6.8771 


AUG-CC-PV5Z 


-9.1289 


-2.2606 


6.8683 


AUG-CC-PCVDZ 


-9.0987 


-2.2371 


6.8616 


AUG-CC-PCVTZ 


-9.1316 


-2.2554 


6.5776 


AUG-CC-PCVQZ 


-9.1293 


-2.2574 


6.8719 


AUG-CC-PCV5Z 


-9.1291 


-2.2607 


6.8684 



of the grid, the grid does not have the shape of the box and there are no basis functions where 
there are no grid points. This means that we are not dealing with box boundary conditions, 
but rather with effective boundary conditions which reflect the shape of the molecule. The 
other point which is not brought out by our explanation is that the BlGDFT grid is adaptive 
in the sense that additional fine grid points are added during the calculation as they are 
needed to maintain and improve numerical precision. 

The implementation of TD-DFT in BlGDFT is described in Ref. P31 . 

8.2. Orbital Energies 

Possibly the most remarkable property of wavelets is how rapidly they converge to the 
basis set limit. Let us illustrate this by comparing highest-occupied molecular orbital 
(HOMO) and lowest-unoccupied molecular orbital (LUMO) energies calculated with DE- 
Mon2k and BlGDFT. The difference of these two energies is the HOMO-LUMO gap, 

^HOMO-LUMO- 

Consider first how deMon2k calculations of A.€homo-LUMO, evolve as the basis 
set is improved (Tabled]) Convergence to the true HOMO-LUMO LDA gap is expected 
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Table 2. Basis set dependence of the HOMO and LUMO energies and of the HOMO- 
LUMO gap (eV) calculated using BlGDFT. 



a Grid spacing of the cartesian grid in atomic units. 
6 Coarse grid multiplier (crmult). 
c Fine grid multiplier (frmult). 

with systematic improvement within the series: (i) double zeta plus valence polariza- 
tion (DZVP) — > triple zeta plus valence polarization (TZVP), (ii) augmented correlation- 
consistent double zeta plus polarization plus diffuse on all atoms (AUG-CC-PCVDZ) — > 
AUG-CC-PCVTZ (triple zeta) -> AUG-CC-PCVQZ (quadruple zeta) -> AUG-CC- PCV5Z 
(quintuple zeta), (iii) augmented correlation-consistent valence double zeta plus polariza- 
tion plus diffuse (AUG-CC-PVDZ) -»• AUG-CC-PVTZ -> AUG-CC-PVQZ -»• AUG-CC- 
PV5Z, (iv) correlation-consistent double zeta plus polarization plus tight core (CC-PCVDZ) 
->■ CC-PCVTZ -)• CC-PCVQZ ->■ CC-PCV5Z, and (v) correlation-consistent valence dou- 
ble zeta plus polarization on all atoms (CC-PVDZ) -> CC-PVTZ -> CC-PVQZ -> CC- 
PV5Z. There is a clear tendency in the correlation-consistent basis sets to tend towards 
values of -9. 13 eV for the HOMO energy, -2.26 eV for the LUMO energy, and 6.87 eV for 
^homo-lumo^ with adequate convergance achieved with the AUG-CC-PVQZ basis set. 

Now let us turn to BigDFT (Table|2]). Calculations were done for several different grids, 
including the high-resolution combination 0.3/8/8 and the low-resolution combination of 
0.4/6/8. Remarkably, except for the very lowest quality grid 0.4/6/8, there is essentially 
no difference between results obtained with the two grids (and even the 0.4/6/8 grid gives 
nearly converged results.) The results are also quite close to, but not identical to those 
obtained with the deMon2k program. The reason for the small differences between the 
converged results obtained with the two programs is more difficult to trace as it might be 
due to the auxiliary basis approximation in deMon2k or to the use of pseudopotentials in 
BlGDFT or perhaps to still other program differences. The important point is that differ- 
ences are remarkably small. 

8.3. Excitation Energies 

Orbital energy differences provide a first estimate for excitation energies. In this case, we 
would expect to see the HOMO — > LUMO excitation at AeuoMO-LUMO ~ 6.9 eV (6.87 
eV for DEMON2K and 6.90 eV for BlGDFT.) A better estimate is provided by the two- 
orbital two-electron model (TOTEM) EH [M EH US for the singlet (S) and triplet (T) 



h g a /m b /n c -(.homo 



(LUMO Af. HOMO — LUMO 



0.4/6/8 -9.0976 

0.4/7/8 -9.1014 

0.4/8/8 -9.1017 

0.4/9/8 -9.1017 

0.4/10/8 -9.1017 

0.3/7/8 -9.1022 

0.3/8/8 -9.1025 



2.1946 6.9029 

2.2028 6.8985 

2.2044 6.8971 

2.2049 6.8967 

2.2049 6.8966 

2.2056 6.8964 

2.2073 6.8950 



34 



Natarajan, Genovese, Casida, and Deutsch 



Table 3. Comparison of lowest excitation energies of CO (in eV) calculated using BlGDFT 
and deMon2k and with experiment. 



State 


BigDFT q 


deMon2k 6 


Experiment 


i 3 sr 


9.84 


9.85 


9.88 


1 3 A 


9.17 


9.21 


9.36 


i J n 


8.94 


8.42 


8.51 


1 3 E+ 


8.94 


8.54 


8.51 


i 3 n 


6.47 


6.05 


6.32 



Q Present work (TD-LDA/TDA) using AUG-CC-PCQZ basis set. 
6 Present work (TD-LDA/TDA) using 0.3/8/8 grid. 
c Taken from Ref. |9"H . 

transition from orbital i to orbital a, 

hwL a = Ae^ a + (ia\fZf-f%/\ai) 

hw^ a = Ae^ a + (ia\2f H + fZf + fZf\ai), (79) 

where 

Ae;_> a = e a - a . (80) 

The TOTEM model often works surprisingly well for small molecules because, unlike 
the Hartree-Fock approximation which is better adapted to describe electron ionization and 
attachment, pure DFT Kohn-Sham orbitals are preprepared to describe excitation energies 
in the sense that the occupied and unoccupied orbitals see the same potential, thus minimiz- 
ing orbital relaxation effects. Inspection of the sizes and signs of the integrals in Eq. d79l ) 
indicates that we should expect, 

hwj^ a < Ae^ a < Huf^ a . (81) 

These is confirmed in Table [3] where the l 3 n and l 1 !! excitations are, respectively, the 
triplet and singlet states corresponding to the HOMO — > LUMO transition. 

Assuming that fx^ a dominates over fxc , we mav even g° a bit further to estimate 
(ia\fn\ai) and (ia\fx^ a \ai) (Fig. El) The calculations are show in Table|4] Comparison of 
(ia\fxc a \ai)^ and (ia\fxi a \ai)^ provides an indication of the quality of the approxima- 
tion of neglecting the {ia\fx£\ai) integral which in this case appears to be excellent. The 
(ia\fn\ai) integrals calculated with the two programs are reasonably close. Interestingly 
the (ia\fx<f*\ai) disagree by about 0.4 eV which, though small, is not negligible. 

Let us now examine the issue of the collapse of the continuum. In Ref. ll92l . it was 
shown that the TD-DFT ionization continuum begins at —euoMO- m exact Kohn-Sham 
DFT, this should be the ionization potential. However typical approximate density func- 
tional underbind electrons and so lead to an artificially-early on-set of the TD-DFT ion- 
ization continuum. This is first illustrated using the deMon2k program and different basis 
sets. Indeed Fig.[8]shows that the states above —enoMO tend to collapse towards —chomo 
rather than converging as they should. This is simply because we are trying to describe a 
continuum which should not be there with a finite basis set. Also seen in the figure is a 
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Figure 7. Estimation of integrals within the TOTEM model. 
Table 4. Estimations of integrals (in eV) within the TOTEM. 
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slight splitting of the l 1 !! excitation energy. This small effect is due to the fact that the grid 
used to calculate xc-integrals in deMon2k has only roughly the symmetry of the molecule. 

Now let us turn to BlGDFT calculations. Figure [9] shows a similar collapse of the con- 
tinuum as the fineness of the grid increases. Interestingly there is no evidence of symmetry 
breaking of the doubly-degenerate l 1 !! state. 



8.4. Oscillator Strengths 

Carbon monoxide is very unusual for small molecules in that absolute oscillator strengths 
have been well studied |93l over a significant energy range and the A l H (l 1 !! in Table [3]) 
is bright and has an accurately determined oscillator strength. See Fig. 2 of Ref. ll93l (as 
well as other references in the same paper) for a graph of measured absolute optical oscilla- 
tor strengths against absorption energy in eV. Table [5]reports our calculated TD-LDA/TDA 
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Figure 8. Singlet and triplet excitation energies for CO calculated using deMon2k 
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Figure 9. Singlet and triplet excitation energies for CO calculated using BlGDFT 



oscillator strengths. As the TDA violates the Thomas-Reiche-Kuhn (TRK) /-sum rule 
RBI it should only be used very cautiously to estimate oscillator strengths. Nevertheless 
the deMon2k value of / = 0.232 is in good agreement with the experimental value of 
/ = 0.1762. As shown in Ref. |[9ll . full TDLDA calculations with asymptotically cor- 
rected potentials give smaller oscillator strengths (0.136 for TD-LDA/LB94 and 0.156 for 
TD-LDA/AC-LDA calculationsEQ.) (Coincidently our own deMon2k full TD-LDA cal- 
cuations without asymptotic corrections give a degeneracyc-weighted oscillator strength of 
0.1752 (bang on the experimental value) but an excitation energy of 8.19 eV.) Since os- 
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Table 5. Comparison of experimental A l Yi energies (eV) and oscillator strengths with TD- 
LDA/TDA experimental A l Ti energies (eV) and degeneracy-weighted oscillator strengths 
(unitless.) 
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0.232 
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a See Table VIII of Ref. lf9TH . 

cillator strengths are quite sensitive to configuration mixing with nearby states, the fact 
that the BlGDFT oscillator strength is larger than the deMon2k oscillator strength may 
be due to the relatively small energy separation between the BlGDFT A 1 II state and the 
artificially-low TD-LDA ionization continuum. 



9. Conclusion 

In this chapter we have tried to give an informative elementary review of a subject largely 
unfamiliar to most theoretical chemists and physicists. Wavelets, once an obscure ripple at 
the exterior of engineering applications, grew to become a regular tsumani in engineering 
circles in the 1990s as the similarity to and superiority over Fourier transform methods for 
multiresolution problems with arbitrary boundary conditions became increasingly recog- 
nized. Though the first applications of wavelet theory to solving the Schrodinger equation 
may be traced back to the mid-1990s ll94l 19311961 . the theory is still not well known among 
quantum mechanicians. Here we have tried to remedy this aberrant situation by trying to 
"make some waves about wavelets for wave functions." 

In particular we have reviewed the theory behind the wavelet code BlGDFT for ground- 
state DFT and our recent implementation of wavelet-based TD-DFT in BlGDFT. Rapid 
progress is being made towards making BlGDFT a high performance computing order- 
N code for applications to large systems. Right now applications to 400 or 500 atoms 
are routine for ground-state calculations with BlGDFT. Our implementation of TD-DFT 
in BlGDFT is by comparison only a rudementary beginning, but it shows that the basic 
method is viable and we are confident that there are no insurmountable obstacles to making 
high performance computing order- N wavelet-based TD-DFT code for large systems. 
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Appendices 

A List of Abbreviations 

For the readers convenience we give a list of the abbreviations used in this chapter in alpha- 
betical order: 

AA Adiabatic approximation. 

APW Augmented plane wave. 

CC Coupled cluster. 

CI Configuration interaction. 

DFT Density-functional theory. 

DM Density matrix. 

FD Finite difference. 

FE Finite element. 

FEM Finite element method. 

GTH-HGH Goedecker-Teter-Hutter/Hartwigsen-Goedecker-Hutter. 

GTO Gaussian-type orbitals. 

H Hartree. 

HF Hartree-Fock. 

HOMO Highest-occupied molecular orbital. 
ISF Interpolating scaling function. 
LCAO Linear combination of atomic orbitals. 
LDA Local density approximation. 
LMTO Linear muffin tin orbital. 
LR Linear-response. 

LR-TD-DFT Linear-response time-dependent density-functional theory. 
LR-TD-LDA Linear-response time-dependent local density approximation. 
LUMO Lowest-unoccupied molecular orbital. 
NS Non-standard. 

MBPT Many-body perturbation theory. 
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MRA Multiresolution analysis. 
S Singlet. 

SOS Sum-over-states. 
STO Slater- type orbital 
T Triplet. 

TD Time-dependent. 

TDA Tamm-Dancoff approximation. 

TD-DFT Time-dependent density-functional theory. 

TD-LDA Time-dependent local density approximation. 

TRK Thomas-Reiche-Kuhn. 

xc Exchange-correlation. 
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